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ABSTRACT 


Finite-element procedures for the Navier-Stokes equations in the primitive 
variable formulation and the vorticity stream-function formulation have been im- 
plemented. 

For both formulations, streamline-upwind/Petrov-Galerkin techniques are 
used for the discretization of the transport equations. 

The main problem associated with the vorticity stream-function formulation 
is the lack of boundary conditions for vorticity at solid surfaces. Here an implicit 
treatment of the vorticity at no-slip boundaries is incorporated in a predictor- 
multicorrector time integration scheme. 

For the primitive variable formulation, mixed finite-element approximations 
are used. A nine-node element and a four-node + bubble element have been 
implemented. The latter is shown to exhibit a checkerboard pressure mode and a 
numerical treatment for this spurious pressure mode is proposed. 

The two methods are compared from the points of view of simulating internal 
and external flows and the possibilities of extensions to three dimensions. 
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CHAPTER 1; INTRODUCTION 


Computational Fluid Dynamics (CFD) is one of the most challenging fields of 
applied sciences. It deals with all the aspects of fiow engineering, from high speed 
reacting fiows around a space vehicle entering the earth’s atmosphere to flows in 
porous media encountered in resevoir engineering. 

The common denominator of this wide range of problems is that their gov- 
erning equations are usually three dimensional, time dependent and nonlinear. 
Moreover, the numerical treatment of practical problems involves the solution of 
very large systems of equations. In recent years, the development of supercom- 
puters combining large in-core storage with very high computing speeds made it 
possible to deal with such problems. Given such hardware capability, the engi- 
neer’s work is the development of numerical procedures adapted to the solution of 
specific classes of problems. 

The present work is a contribution to the development of numerical schemes 
for the solution of the Navier-Stokes equations. 

The use of the vorticity stream-function form of the Navier-Stokes equations is 
appropriate for the study of two-dimensional flows [1,2]. For two-dimensional sit- 
uations, the vorticity stream-function equations reduce to a set of scalar equations 
so that the number of degrees of freedom of the problem is only two. Another 
main interest of that formulation is that the incompressibility constraint is au- 
tomatically satisfied by the computed fiow field. However, solving the vorticity 
transport equation presents difficulties due to the lack of boundary conditions at 
solid surfaces. A feature of this study is the proper numerical treatment of the 
boundary conditions for the vorticity at solid surfaces. A mixed finite-element for- 
mulation is used that allows the implicit computation of the trace of the vorticity 
on no-slip boundaries at each time step. 

The primitive variable formulation of the Navier-Stokes equations is suited 
for two-dimensional and three-dimensional calculations, and the treatment of the 
boundary conditions presents no difficulties [3,4]. The major problem associated 
to that formulation is the necessity of using different orders of interpolation for the 
velocity and the pressure [5,6]. In the present study mixed finite-element interpo- 
lations using quadratic approximation for the velocity and linear approximation 
for the pressure have been used. A four-node + bubble element, and the nine-node 
Lagrange element have been tested. 

A common problem associated to both formulations is the treatment of the 
advection term . It is well known that the Galerkin approximations, which lead to 
centered discretization of the convection term, may create numerical oscillations 
for convection-dominated problems especially in the presence of sharp gradients in 
the solution [7]. On the other hand, such schemes introduce no numerical dissipa- 
tion. The schemes introducing streamline diffusion via a suitable Petrov-Galerkin 
formulation are more robust than the centered schemes; they behave nicely at 
sharp fronts but of course involve relatively more numerical dissipation especially 
for high wave numbers. Such procedures have been applied to various fluid dy- 
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namics and convection-difTusion-reaction problems [8-12]. In the present study, we 
employ streamline-upwind/Petrov-Galerkin (SUPG) formulations for the solution 
of the vorticity transport equation in the vorticity stream-function formulation and 
for the solution of the momentum equation in the primitive variable formulation. 
The weighting functions are based on temporal as well as spatial discretizations 
( 11 ). 

The present study has been devoted to the development of software for two- 
dimensional simulations. However, in the case of the primitive variable formu- 
lation, the selection of procedures and the design of the code are such that the 
extension to three dimensions would present no major difficulties. In the case of 
the vorticity stream-function formulation, the simplicity of the equations is re- 
lated to the dimensionality of the problem and no extension to three dimensions 
is foreseen. 
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CHAPTER 2; PROBLEM STATEMENT 


2.1 The velocity-pressure formulation 


Let n be a domain of and T be a positive real number. The flow of a viscous 
incompressible fluid in fl x [0,T] is described by the Navier-Stokes equations: 
Momentum equation 

/fll 

— + (u.V)u = V.O in fl X [0,T] (2.1) 

ot 

Continuity equation (Incompressibility constraint) 

V.u =0 on n (2.2) 

where u(x,t) = {«i(x,t)}?^i is the velocity fleld, o(x,t) = is the 

Cauchy stress tensor. For simplicity, we assume that the density of external forces 
is zero. For a Newtonian fluid of unit density, the stress tensor is given as 


a = -pi + 2i^e(u), (2.3) 

where p is the pressure and r(u) is the tensor of rates of deformation, i.e, 

ff(«) = + (^«)^) (2-4) 

We assume that the boundary P of fl is sufficiently smooth and admits the 
following decomposition : 

T=TgUTH (2.5) 

and 

r,nPfc = 0, (2.6) 

where Pj and P^ represent the Dirichlet- and Neumann- type boundaries. Hence, 

u(x,t) = g(x,t) on Pg X [0,r] (2.7) 

n.cT = h(x,t) on Th x [0,T], (2.8) 

where n is the outward unit normal vector to P/,. 

In the case of a pure Dirichlet problem ( Pa = 0), the specified boundary 
velocity must satisfy 

j g.n <fP = 0, (2.9) 

i.e, the net flux of matter through the boundary should be zero. 

The initial condition associated to this problem is 


u(x,0)=Uo(x) on n 


( 2 . 10 ) 
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For the we!l>posedness of the problem, the initial field must be divergence-free, 
i.e, 

V.uo = 0 (2.11) 

We will see later the importance of the discrete counterpart of that condition in 
the solution process. 

The system of eqns (2.1)-(2.2) and the associated boundary and initial con- 
ditions are called the velocity-pressure (or primitive variable) form of the Navier- 
Stokes equations. 

2.2 The vorticity stream-function formulation 


The vorticity transport equation can be obtained by taking the curl of the 


momentum equation, 

V X {^ + (u.V)u} = V X {- Vp + i/V*u} (2.12) 

Now, we use the vector identity [13], 

(u.V)u = ^V|u|* - u X (V X u) (2.13) 

Replacing eqn (2.13) in eqn (2.12) and using the relations 

V X V|u|*= 0 (2.14) 

V X Vp = 0 (2.15) 

V*u = -V X (V X u) (2.16) 

( in eqn (2.16), we have used V.u =0 ), we obtain 


d{V X u) 
dt 


-Vx(uxVxu) = -j/V X V X V X u 


(2.17) 


The vorticity is defined as 

w = V X u (2.18) 

By substitution of eqn (2.18) in eqn (2.17), we obtain the vorticity transport 
equation, 

~ + (u.V)u; = (w.V)u -I- (2.19) 

ot 

In the case of a two-dimensional flow, eqn (2.19) simplifies to 


— -|- u.Vu? = uV^u} 


( 2 . 20 ) 
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The velocity u being solenoidal, there exists a vector field 

® = (0,0, '9) such that 

u = V X ® 

(2.21) 

The nonzero component of 9 is called the stream function. A relation between 

the vorticity and the stream function is obtained as 


u; = V X Vx® 


= 

(2.22) 

The vorticity stream-function form of the Navier-Stokes equations for two- 
dimensional fiows is the set of coupled scalar equations given below : 

dijj o 

+ u.Vw = i/V^w on n X [0, T] 

(2.23) 

= onnx[0,r] 

(2.24) 

where n is a domain of R*, u(x,t) is the velocity, w(x,t) 

is the vorticity, ^(x,t) 

is the stream function, and u is the kinematic viscosity, 
stream function are related to the velocity field through 

The vorticity and the 

duj dui 

(2.25) 

dxi dx2 


(2.26) 

“ dX2 


(2.27) 

U2 = - — 
dxi 

We assume that the boundary F of the domain H has the following decompo- 

sitions: 


r = Fj U Fg U Ffc 

(2.28) 

0 = Fj n Ffc 

(2.29) 

0 = FcnFfc 

(2.30) 

0 = F, n Fg 

(2.31) 

F = Fj U Fg U Fjr 

(2.32) 

« = rjnr» 

(2.33) 

0 = FGnF^ 

(2.34) 

0 

c 

II 

(2.35) 
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where and Th represent Dirichlet- and Neuman-type boundaries for 'i', Fj 
and F^ represent Dirichlet- and Neumam-type boundaries for u;, and Fc represents 
both Dirichlet- and Neuman- type boundaries for ^ . Hence, 


«'(x,t) = '*',(x,t) 

onF, X [0,r] 

(2.36) 

«'(x,0 = ^g(x,0 

on Fg X [0,r] 

(2.37) 

n.V^'(x,t) = Ur 

on Fg X [0,r] 

(2.38) 

n.V’i'(x,t) = h{x,t) 

onFfc X (0,r) 

(2.39) 

w(x,«) =ug(x,t) 

on Fj X [0,r] 

(2.40) 

un.Vu(x,t) = h(x,t) 

on Fji X [0,r] 

(2.41) 

The boundary F^ represents the solid surfaces in the flow domain, 
condition for vorticity in eqn (2.23) is given as 

The initial 

w(x,0) =wo(x) 

on n 

(2.42) 


Remarks : 

1) In the case of a two-dimensional flow, the vorticity-stretching term (w. V)u 
is zero and the vorticity transport equation contains only one nonlinear term. This, 
together with the scalar nature of the equations, is the reason why the vorticity 
stream-function formulation is very convenient in two dimensions. 

2) The vorticity transport eqn (2.23) belongs to a large class of convection- 
diffusion problems . In this case, the transported quantity is the vorticity or. 

3) The flow field obtained from the solution of eqns (2.23) and (2.24) is 
divergence-free by construction. 

4) For flow domains having solid boundaries ( e.g., walls or obstacles ) the 
no-slip and no-penetration conditions yield constraints on the derivatives of 'I', i.e, 

n.V'i' = Ur (2.43) 

f.V'i' = 0 (2.44) 

where r and n are the tangential and normal unit vectors at the surface whereas Ur 
is the tangential component of the velocity at the wall. Moreover, the stream 
function assumes constant values along solid boundaries. However, it is difficult to 
derive any boundary condition for vorticity at solid surfaces. The boundary values 
of the vorticity at solid surfaces are related to the local boundary-layer profiles 
and are thus time dependent. On the other hand, they do not arise naturally 
from the no-slip and no-penetration conditions. A classical procedure to update 
the wall values of u in multiply-connected domains is to include the solution of 
an auxiliary Poisson problem in a semi-implicit time integration scheme for eqns 
(2.23)-(2.24) [l]. An alternative approach to the numerical treatment of vorticity 
at no siip-boudaries is discussed in [2]. 

For flows in symmetrical domains having only one inner boundary, we propose 
[14] to take advantage of the time discretization of the problem to treat the wall 
values of u; as unknowns in the discrete Poisson equation arising from eqn (2.24). 
This process will be described in Chapter 3. 
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CHAPTER 3: NUMERICAL PROCEDURES FOR THE 
VORTICITY STREAM-FUNCTION FORMULATION 

3.1 Variational formulation 

In the sequel, we will use the Hilbert space 

^ 2 ( 0 ) = {$1$ is possibly discontinuous and / |$|*dn < + 00 } (3.1) 

Jo 

and the Sobolev space 

iTi(n) = {«|« € Lj(0), ^ 6 1’(n),v.- = 1,2}. ( 3 . 2 ) 

We introduce the following spaces: 


^{n) = {w|iy € Hi{tl),w = 0 on F^} 


(3.3) 


S(n) = € Hi{n),« = on on Fc} 


(3.4) 


V{n) = € Hi(n),u; = 0 on Fj U Tq} 


(3.5) 


S(n) = {w|u; € lfi(n),w = Uj on Fj,w = ug onFo) 


(3.6) 


The variational form of the vorticity stream-function problem is the following: 
Find u in <S(Q) such that for all w in V, 

f u;(^ +u.Vu;)(in+ f uVw.VudU= f whdT (3.7) 

Jo Jo Jrf^ 

and find 'if in S(0) such that for all w in V(n), 

I Vw.V'ifdVl — I wudil = j whdT + [ wurdT (3.8) 

Jo Jo JTh Jva 
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3.2 Finite-element formulation 

Let denote a finite-element discretization of the computational domain fl 
into subdomains 0* , e = 1, 2, ne/, where nel is the number of elements. Let M 
be a vertex of Q* . We assume 

nel 

(3.9) 

e=l 

nel 

0 = Pi n* ( 3 . 10 ) 

e=sl 

We associate to the finite dimensional spaces, 

^f(n) = € c'®(n),$*in« € p‘,vn* € (3.11) 

where is the space of first order polynomials in xi,X 2 , 


V'^(n) = = 0 on T,}, 


(3.12) 


5 *(n) = € /ff (n), on Tg, c- on r©}, ( 3 . 13 ) 


v'‘(n) = € pf (n), 0 on Tj U To), (3.14) 


5^(n) = {w*|u;* € ff*(n),u;* Ug on TyjO/* cs wg onFc}. (3.15) 

In the above definitions, we assume that 9g, ufg and ufc ere known quantities. 
The discrete variational problem associated with eqn (3.15) is given as follows: 
Find w* in 5*(n), such that for all w* in 



u.Vu;^)dn + f 

in 


i/Vw^.Vuj^dn+ 


+ u.Vw* — i/V 


^u>^)dn= f w^hdT 

Jtl 


(3.16) 


where 5 is a C“^(fl) perturbation to the weighting function w^. The selection of 
the Petrov-Galerkin perturbation will be further addressed in Section 3.4. Based 
on u;* ( solution of eqn (3.16)) we define the following space: 

5^(n) = € P(‘(n),u;^(M) = u,^{M) if M ^ To) (3.17) 
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The discrete variational form associated with eqn (3.8) can be stated as follows: 
Find ’i'* € 5^, and w* € Sq such that for all tw* 6 

f - f = f w^hdT + f w^UrdT (3.18) 

Jci J(i JTk Jra 

Remarks : 

1) Note that if 6 = 0 we have a Galerkin formulation; if £ 0 we have a 

Petrov-Galerkin formulation. The perturbation weights only the element interiors 
and does not affect the boundary terms. 

2) Brooks and Hughes [8] have shown that in the case of general convection- 
diffusion problems, for rectangular elements , the perturbation does not affect the 
weighting of the diffusion term, and for reasonable element shapes the contribution 
is expected to be negligible. 


3.3 Temporal discretization 

The system of ordinary differential equations obtained from eqns(3.16)-(3.18) 
is 

Ma + Cv = F (3.19) 

v(0) = Vo (3.20) 

Kd - Mv* = E (3.21) 

where v is the vector of nodal values of and a is its time derivative, whereas d 
is the vector of nodal values of and v* is the vector of nodal values of In 
eqns(3.19)-(3.21), M is the mass matrix, C is the sum of the diffusion matrix and 
the linearized convection operator, K is the discrete Laplace operator ; F and E 
are generalized force vectors containing the fluxes and the specified values at the 
boundary. 


This system is solved with the following semi-implicit predictor /multi- 
corrector algorithm [8]: 


i) Initialization: Vn,an«dn are known. If n = 

0, the initial conditions are 

used. 


ii) Prediction: 


= Vn-h(l-'7)Atan 

(3.22) 

a^+1 = 0 

(3.23) 

Kd^^i — = En+i 

(3.24) 

„o _ 

dz, ’ azi > 

(3.25) 
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Here At is the time step and 7 (0. < 7 < 1.) is a parameter governing the stability 
and accuracy of the integration scheme. 


iii) Correction: 


M Aa„^i = Fn+i — MaJj^.1 — 

(3.26) 

with M* = (M + 'fAtC) 

(3.27) 

KVi =<‘1 +^^<^4+1 

(3.28) 

“ ®n+l + ^®n+l 

(3.29) 

Kdii*. -Mv‘+i. =E,+. 

(3.30) 

di, ’ dx, ’ 

(3.31) 


The iterations continue until the norm of the residual becomes smaller 

than a predetermined error tolerance. 

3.4 Selection of the Petrov- Galerkin perturbation 

In the present study, we have used a streamline-upwind/Petrov-Galerkin pro- 
cedure. In this formulation, the perturbation to the shape function is given as 

6 = CirZ^S.VNa, (3.32) 


where 


C2r is the algorithmic Courant number [11], 

u 

8 = 7^ — T is a vector of unit length pointing in the advection direction, 

il«ll 

h is the element length, 

s is a function of the element Peclet number Pe given as [8] 


z 



1 


if 0 < Pe < 3 
if Pe > 3 


where Pe = ||u|l— is the element Peclet number. 

Various selections for Ctr and their stability /accuracy properties are discussed 
in [llj. In the present study, we have used 


Cir = 1 


(3.33) 


and 


C2r = CAt 


(3.34) 
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where 


(3.35) 


is the element Courant number. 

These selections of the Petrov-Galerkin perturbation ore such that the scheme 
introduces a certain amount of artificial diffusion in the convection direction, but 
induces minimal such diffusion in the cross-convection direction. This method 
can be recast in terms of the use of an additional anisotropic diffusion tensor [8] 
added to the physical one which is isotropic in the case of Newtonian fluids. Such 
procedures have also been introduced through the use of a generalized governing 
equation [15] obtained by subtracting from the original differential equation the 
scalar product of its gradient by a vector of free parameters associated to each of 
the coordinate directions. 

All the so-called Petrov-Galerkin formulations are interided to prevent the oc- 
curence of oscillations in the computed flow field. Such oscillations in the solution 
can be obtained when regular Galerkin formulations are used. Typically, in the 
presence of shocks or boudary layers, if the mesh is not fine enough to resolve 
the discontinuity and sharp gradients, the solution might be oscillatory. Usually, 
the oscillations in the solutions can be avoided by a mesh refinment where the 
mesh spacing is set smaller than the boundary layer thickness [16] . However, this 
approach can lead to a dramatic increase in the number of nodal points. The size 
of the problem is then likely to become very large, especially for three-dimensional 
simulations. The use of Petrov-Galerkin formulations is an alternative to mesh 
refinment for obtaining accurate wiggle-free solutions in a cost-effective way. 
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CHAPTER 4: NUMERICAL PROCEDURES FOR THE 
PRIMITIVE VARIABLE FORMULATION 

4.1 Variational formulation 

Using the definition of H\ (H) introduced in Section 3.1, we define the following 
spaces: 


v(n) = {w|w 6 ^i(n)xB-i(n),w = o on r,} (4.1) 

5(n) = {u|u € ^i(n) X jBTi(n),u = g on r,} (4.2) 

The variational form of the problem defined by eqns (2.1)-(2.2) can be stated 
as follows : Find u € S{tl) and p € such that for all w € V(n) and for 

allgeHi(n), 


f w^dn + f w(u.v)udn = [ wv.odn 
Ja ^ Jn Ja 


(4.3) 


and 


I gV.udn = 0 
Ja 

Note that 

I wV.ffdCl = / V.(w<T)dn — / <T.Vwdn 

Ja Jn Ja 

= / w.(n.o)dT — / ff.Vwdn 

Jr Ja 

So eqn (4.3) can be rewritten as 

I w^dn+ / w(u.V)udn = f w.(n.< 7 )dr - / tr.Vwdn 
Ja ^ Ja Jth Ja 

Moreover, 

a.Vw = — pI.Vw + 2i/c(u)Vw 
= — pV.w + 2t/£(u)e(w) 


We obtain 


f w^dn + f w(u.v)udn — f pV.wdn 
Ja ^ Ja Jn 

+ / 2 t/c(w)c(u)dn = / w.(n.cr)dr 

Jn Jtk 


(4.4) 


(4.5) 


(4.6) 


(4.7) 


(4.8) 
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Using the boundary condition (2.8), we get 


f w^dn + f w(u.V)udn — f pV.wdn 

Ja ^ Jn Jn 

+ / 2t/c(w)c(u)dn = / w.hdr 
Jn Jtk 


(4.9) 


The variational form of the continuous problem (2.1)-(2.2) is given by eqns 
(4.4) and (4.9). 

The initial condition associated to the variational problem is 


/ w{u(x,0) - Uo(x)}dn = 0 

Ja 


(4.10) 


4.2 Finite-element formulation 

Let denote a finite-element discretization of the computational domain fl 
into subdomains O*, e = 1, ..., nel where net is the number of elements. We assume 


nel 

(4.11) 

e=l 


nel 


0= n 

e=rl 

(4.12) 


To obtain a discrete variational problem, the functions w,u,p, q are approximated 
by piecewise polynomial functions so, we introduce the following 

finite dimensional spaces: 

V^^(n) = {w'^lw* € ffi(n) X ffi(n),w'‘|n« € = 0} (4.13) 

= {u'^lu'^ € ^Ti(n) X fTi(n),u*|n« € P’*,u*|r. = g'‘> (4.14) 

ffim = {9*1/ € fri(n),/|n. € (4.15) 


where P“ (resp. P^) is the set of piecewise polynomial functions of order a (resp. 
fi). The selection of the order of the polynomial interpolants will be discussed in 
the sequel. 


13 


The Petrov-Galerkin variational form of eqns (4.9) and (4.4) is the following: 
Find € S^{a) and € ^{‘(n), such that for all w^‘ € ^'‘(n) and for all 

Z 


[ [ w*(ii^.Vu*)<(n + 2i/ / ff(w^)ff(u*)(in — [ p*v.w*<in+ 

Jn ot Jq Jq Jq 



+ u.Vu^ — Vcr*)dn = / w^h<fT 

Jvk 


f q^V.n^dn = 0 
Ja 


(4.16) 

(4.17) 


where ^ is a C"^(n) perturbation. The selection of the Petrov- Galerkin pertur- 
bation will be discussed later. Note that if ^ = 0, we have a Galerkin formulation; 
if 5^0, we have a Petrov-Galerkin formulation. 

We can rewrite eqn (4.16) as 


f [ w^u^.Vu^dfl + 2i/ / e(w*)c(u*)<in 

Jn ^ Jn Jn 

-L 


p^v.w’^dn 


“£/ = J w'‘hdr 


(4.18) 


where + 5. In the sequel we will assume that the contribution of the 

Petrov-Galerkin perturbation to the stress-divergence term ( the last term on the 
left-hand-side) is negligible. 


4.3 Temporal discretization 


4.3.1 Time-integration scheme 

The discrete eqns (4.16)-(4.17) have the following equivalent matrix form: 

Mv + N(v) + Cv - Gp = F (4.19) 

G^v = D (4.20) 

where v is the vector of unknown nodal values of the discrete velocity field, v is 
its time derivative and p is the vector of unknown nodal values for the discrete 
pressure field. M is the mass matrix, N(v) is the nonlinear convection vector, C is 
the diffusion matrix, G is the discrete gradient operator. F and D are generalized 
force vectors containing the prescribed fluxes and displacements. 

The time integration of this system is obtained through the following 
predictor-multicorrector algorithm {8|: 
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1) Initialization: 


V«, ❖n, Pn are known. 

If n = 0, the initial conditions Vo and Vq are used. In the continuous problem, 
there is no initial condition required for the pressure, and neither is there any in 
the discrete problem. However, as will be shown later, the initial discrete velocity 
field mtist satisfy the approximate incompressibility constraint. 


G^vo = 0 

(4.21) 

2) Prediction : 


V°+l = Vn + (1 - 7) AtVn 

(4.22) 

= 0 

(4.23) 

where the definitions of At and 7 are the same as in Section 3.3. 
3) Correction : 


+ N(vi;>,) + cvjy, - = f„+, 

(4.24) 

= Di*+1 

(4.25) 

We introduce the incremental forms given below: 


= '^n+l + 

(4.26) 

vi+\ = v‘n+l + ^<+1 

(4.27) 

PiVl = Pn+1 + ^Pn+1 

(4.28) 


Substituting eqns (4.26), (4.27) and (4.28) into eqn (4.24) yields 


+ N(v ~ ^Pn+1 

+ MAvi+i + CAvi+i - = Fn+i (4.29) 

A Taylor expansion of the vector N(v*,^.i + Av|,^i) in the vicinity of 
yields 

N(v;+i + Av*„+i) c- N(v;+i) + ^ (4.30) 

From eqns (4.29) and (4.30) we obtain 

MAv*„+i + + CAv*„+i - GAp*„+i 

= Pn+i - + n-yUi) + - GpUi) (4.31) 
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Let 


R*»+. = F»+i - {M^Ui + N(vU.) + Cvi+i - Gp’„+.} 

On the other hand, 

=iA(A<+,; 

SO, eqn (4.31) can be rewritten 

[M + + C)] - GApU, = RUl 


The discrete continuity equation reads 


G^Av‘ 1 = Dn+i - G V 


n+l 


or 


Let 


Let 




Bi+l = - G'vUi) 


A = M + iA1(^L +C) 


dv 'n+i 

The residual form of the equation system (4.19)-(4.20) is 


(4.32) 

(4.33) 

(4.34) 

(4.35) 

(4.36) 

(4.37) 

(4.38) 

(4.39) 


G’^ AvUi = Bi+i (4.40) 

If AvJj+i and ApJ^^i are obtained by solving eqns (4.39)-(4.40), then the scheme 
is implicit. In the present formulation, an explicit scheme is used for the discrete 
momentum equation. Let 

M* = (4.41) 

be the lumped counterpart of M obtained by a column-sum technique. The ma- 
trix M* replaces A in eqn (4.39) and we end up with the system, 

M‘A<+, - G'f ApU, = RUi (■<■«) 

G’'A*U, = BUi (4.«) 

The following change of variable is introduced [8], 

A*;V: = - M-'GApi+, (4.44) 
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Eqn (4.42) reduces to 


(4.45) 


M*Av;Vi = RUi 

Replacing eqn (4.44) in eqn (4.43), we obtain an equation for as 

= Bi+i - G’A^;V, (4.46) 

Now, let K = G^M* *G. The pressure increment is obtained from 

KAp‘„+, = - G»Ay;‘+, (4.47) 

The complete acceleration increment is recovered from eqn (4.44) through 

A<h.i = + M-'GAp‘„, (4.48) 

The discrete velocity and pressure fields are obtained from eqns (4.27)- (4.28). 

The use of mass-lumping has many advantages. The left-hand-side matrix 
becomes diagonal and the storage required is reduced to a vector of length equal 
to the number of equations. Moreover the solution of the momentum eqn (4.45) 
can be obtained through a scalar product of vectors. This is an important asset 
in the framework of vectorization, because vector manipulations are among the 
fastest operations on supercomputers. Of course the use of an explicit scheme is 
expected to reduce the convergence rate of the solution. However, the increased 
number of iterations is fully compensated by the gain in storage. On the other 
hand, the matrix K is symmetric and the solution of the pressure equation can 
then be obtained through conjugate gradient type techniques. Note that K is not 
obtained by a classical assembly of element-level matrices but rather by mapping 
of element-level products into their location in the global matrix. 

The reduction in storage requirements introduced by the use of an explicit 
treatment of the momentum equation is very important regarding the extension of 
the code to three dimensions. In that case, the study of practical problems involves 
solution of systems too large to allow the use of implicit algorithms. Owing to 
the implicit nature of the pressure, the discrete continuity equation cannot be 
solved by an explicit technique. However, the storage required for the matrix K is 
small with respect to that required for A ( if the momentum equation was solved 
implicitly ) because the pressure is a scalar quantity and because the order of 
interpolation of the pressure is less than that of the velocity. There are far fewer 
unknowns for the pressure than for the velocity and the size of K is much smaller 
than the size of A. 

4.3.2 Selection of the initial condition 

From a mathematical point of view, for the well-posedness of the problem, it 
is fundamental to have a divergence-free initial flow field [5], ie, 

V.uo = 0 (4.49) 
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From a numerical point of view, this condition translates to 

G^vo = 0 (4.50) 

that is the discrete initial flow field should satisfy the discrete incompressibility 
condition. Failure to satisfy this constraint can lead to nonconverging solutions. 
In fact, suppose 

G^Vq ^ 0 and Vo = 0; 
then, the prediction step gives 

V® = Vo (4.52) 

V? = 0 (4.53) 

After the first iteration, the discrete velocity field is 

v} = V? + + M-^GAp}) (4.54) 

where Av^^ is the solution of the discrete momentum equation at the first itera- 
tion, and Ap} is the unknown pressure increment. Because we require G*^vJ = 0, 
we obtain 

GTvJ = G^'v? -h 7AtG'*‘Av;‘ + ^AtG^^M'^^GApJ = 0 (4.55) 

or 

^0 


-1 


0 A • 

, - G^Av, 


= — ^G^vo-G^AV 


•1 


(4.56) 


So if G^Vo ^ 0, the right-hand side of the discrete Poisson equation varies like 
At~^. This can produce the nonconvergence of the algorithm. However, it is 
not easy to specify an initial flow field satisfying the discrete incompressibility 
constraint. In the present study, we started the computations with an arbitrary 
convenient initial condition, typically Vo = 0 except on the Dirichlet boundaries. 
After one time step, the solution is such that G^Vi = 0, but the pressure field 
is ususally erroneous. So because we have no initial condition for p, we set it to 
zero and use Vi as an effective initial condition. 


4.4 Selection of the Petrov- Galerkin perturbation 

As for the vorticity stream-function formulation, we have used a streamline- 
upwind/Petrov-Galerkin scheme. In this scheme, the perturbation is again of the 
form, 

6 = CirZ^a.VNa, (4.57) 

with the definitions of Section 3.4 except for z. Here, we have taken 

s = coth(— ) - — (4.58) 

This expression is an approximation of the formula given in [15] for the nine-node 
Lagrange element in two dimensions. 
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4.5 Selection of a mixed finite-element interpolation 


4.5.1 Problem statement 

In the primitive variable formulation of the Navier-Stbkes equations, the or- 
ders of interpolation for velocity and pressure a and 0 cannot be selected arbi- 
trarily. The problem can be set in the following terms: 

The discrete equation, 

Gp = 0, (4.59) 

should have the unique solution p = 0. Here p is the vector of nodal pressures, 
and G is the discrete gradient operator. Recall that the construction of the discrete 
gradient operator involves products of functions coming from the pressure and 
velocity spaces. The nonzero solutions of eqn(4.59), if they exist are called spurious 
pressure modes. They can be picked-up in the solution process and pollute the 
pressure field. However their presence does not affect the velocity solution. In 
fact, assume that ( u,p) is one solution of the discrete system (4.19)- (4.20), then 
( u, p + p«) with p« ^ 0 satisfying 6p« = 0 is also a solution. 

Equal orders of interpolation for velocity and pressure, i.e, a = 0 are known 
to produce such spurious pressure modes. The usual cure is to use mixed finite 
-element interpolations, typically by interpolating the velocity with a polynomial 
function of degree higher than that used to interpolate the pressure. The two 
most usual combinations are bilinear velocity ( a = 1) /constant pressure ( /d = 0) 
and quadratic velocity ( ot = 2) /bilinear pressure ( )? = 1). However, the bilinear 
velocity /constant pressure formulation has been used with mixed success [6,17]. In 
particular, it has been reported to produce spurious pressure modes under certain 
boundary conditions. 

For the present work, we have chosen to use a quadratic velocity /bilinear 
pressure formulation. 

4.5.2 The four-node + bubble element 

The first element we have used in the computations is a bilinear + bubble 
velocity /bilinear pressure element. It contains five nodal points for velocity, one at 
each of the four vertices of the quadrilateral and one at the center of the element. 
The four nodes for the pressure are at the four vertices of the element (Figure 
4.1). This element is not usual for the solution of the Navier-Stokes equations. 
Our goal when selecting it was to use a mixed interpolation of the type quadratic 
velocity/bilinear pressure while reducing the cost associated with the computation 
of the element level arrays. In fact, in that situation, the number of velocity 
unknowns is 10 for each interior element; thus, the element level arrays such as 
mass, diffusion and convection matrices are 10 x 10 matrices. The element level 
gradient operator is a 10 x 4 array. However, this element is only an extension of a 
bilinear velocity/bilinear pressure element and is not totally free of the drawbacks 
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master element (b) 


2 ( 




associated with the equal order of interpolation. 

The shape functions in the parent domain are (with the numbering of Figure 
4.1) given below: 

For velocity, 


(4.60) 

JVa = j(l + «)(!-»)- 

(4.61) 


(4.62) 


(4.63) 

jv. = (1 - e’)(i - ,») 

(4.64) 

The interpolation of velocity is done through an incomplete second-order polyno- 
mial. 

For pressure, 

1 

1 

wH 

r-i 1 

II 

(4.63) 

II 

1 

+ 

1 

(4.66) 

J»3 = j(l + {)(! + >;) 

(4.67) 

i». = i(l - {)(! + „) 

(4.68) 


4.5.3 Existence of a checkerboard mode for the bubble*element 

We will show that a checkerboard mode can exist when the mixed interpola- 
tion uses the bubble-element. 

Consider a regular mesh of rectangular elements (Figure (4.2)). Using the 
local numbering of the nodes of element e, the pressure shape functions expressed 


in terms of the physical coordinates are 

Ni{x,y) = ^(i,+i - x)(y,+i - y) (4.69) 

■^^(ic.y) = - ®)(y - yy) ( 4 - 72 ) 
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X 


Figure 4.2. Global node numbering for a regular mesh. 
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The correspondance between the global coordinates (z, y) of a point M in 
the physical domain and the coordinates (^,^ 7 ) of its image obtained by mapping 
the physical element into the parent (or master) element is given by 




The element level discrete form of a checkerboard can be expressed as 

p: = p*(Ma) = k(-l)‘ (4.75) 

where Ma is a vertex of element e. Without loss of generality, we take A: = 1 . 
Using the expansion of p‘ over the shape functions associated to element e, the 
interpolation of the spurious pressure mode is 

p‘ = E 

a=l 

= -iVf + -Ni + Nl (4.76) 






Xj + Xj+I ^ 

2 2 ^ 




aasl 


yj + iij+i * 


(4.73) 


(4.74) 


Using the above relations, this becomes 


(4.77) 


According to eqn (4.59), a discrete spurious pressure p mode can exist if it 
satisfies 

Gp = 0 

The entry Ai of vector Gp is given by 

net . nenp 

{Gp)x, = E/ ^•.<(E^‘'’«‘'« 

e*l*'f*' 6=1 

nel f + l ^+1 

= E y_, + JV.,,„,i)(-£n)j<i«<i-» ( 4 . 78 ) 
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where a is the local node number corresponding to A on element e, nenp is the 
number of element nodal unknowns for the pressure, ne/ is the number of elements 
and j is the Jacobian determinant of the mapping from the physical domain into 
the parent domain. 


We can express and iVa,,, as follows: 

For a = 1,2, 3, 4, 


= Jca(l + »?a»7) + ^^(l ~ »7*) 

(4.79) 

■/V^a.n = ^^a(l + ^aO + ~ 

(4.80) 

For a = 5, 


JV.,< = -2f(l - 8’) 

(4.81) 

JVs,„ = -2ll(l - {’) 

(4.82) 


By substitution of eqns (4.79)-(4.80) into eqn (4.78), we obtain, 
for a = 1,2, 3, 4, 


^ r+i r+i 1 1 

{Gp}x, = ^ y ^ j ^ (-|o(l + -^(1 - 

^+1 C+l 1 1 


11 11 

= X^y ^ y ^ ~ + -^'^ri^)U3didri 

r+1 r+i 11 11 

S y J y J - -Una^^rj - 2^»?^ + ^fn^)n,ijd^drj 

for a = 5, 


(4.83) 


««* .+1 «-+i 

{Gp}a. = 5^ / / (-2^(1 - ri^){-(T})Kdd^dT} 

^+l r + l 

+ Xj y J y J (“2r/(l - t^){-^v)v,ijd^dri 
««< /•+! /•+1 

= ^y ^ y ^ 2(^*q - (^rj^)^,ijd^dT} 

/•+! /•+! 

+ y ^ y ^ 2 (^» 7 ^ - (^r)^)n,ijdCdv 


(4.84) 
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For a regular mesh of rectangular elements, and j are constant. Thus 

all the terms under the integral signs in eqns (4.83) and (4.84) contain odd powers 
of either ^ or rj. Their integrals on the domain centered on (0,0) are equal to 
zero. This means that the checkerboard pressure mode satisfies the equation 


Gp = 0 


4.5.4 Treatment of the spurious pressure mode 

The correct pressure field can be recovered through a filtering technique [I8j. 
In the presence of a checkerboard pressure mode, the computed pressure field can 
be written 

P = P + Pc (4.85) 

where p is the correct pressure field and pc is the checkerboard pressure. For a 
regular mesh of rectangular element, the interpolation of the checkerboard pressure 
over one element is 

p\{Ma) = (-l)“fc (4.86) 

where Ma is a vertex of the pressure discretization, and k is an arbitrary constant. 

To filter the pressure solution and recover an acceptable pressure field, we use 
the following two-step approach: 

First step: Averaging. 

The pressure p is averaged over each element through 

?'(«■) = 

a=l 

where Pa is the value of the pressure at node a of element «. In the case of a 
regular mesh of rectangular elements, we have 

, * 

^ (4.88) 

^a=l 

that is , the averaging removes the checkerboard. 

Second step: Least-squares approximation. 

We recover the nodal values of the pressure through a least-squares approach. 
Let p denote the unknown discrete pressure field, and let p denote the discrete 
averaged pressure field. Let J(p) be the functional 

•^(p) = / (P - (4.89) 

Jn 
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Let us expand p over the basis of the pressure approximation as 



np 

p = ^ PbNb 

(4.90) 


B=l 


where np 

is the number of unknown nodzd pressure values. Then J(p) 

can be 

written 

. np 



^{p) = 1 ip- PB^B^dii 

(4.91) 




We look for p such that J(p) is a minimum, i.e., 



o 

II 

(4.92) 

But 

dJ 

6J = 0=^ — =0 VC = l,...,np 

opc 

(4.93) 

Now 

dJ f 

5 — = -2 / (p- y; vbNb)Nc«i 
opc Ja 

(4.94) 

and 



dJ 

dpc 

= 0 =► / pNcdU = / (i^ PBNB)Ncdn VC = l,...,np 

•'o Jn 

(4.9S) 

or equivalently 



np . . 

(5Z / NBNcdn)pB= / pNcdil VC = l,...np 

•'n •'n 

(4.96) 

Let p be the vector of unknown nodal values of p. Then the above equation 

system can be written as 



Mp = F 

(4.97) 

where 

Mcb = { / NBNcdn} 
Ja 

f 

(4.98) 

and 



Fc = / pNcdtl 
Jo 

(4.99) 


The matrix M has the structure of a mass matrix arising from a finite element 
discretization. 


Remark: 


26 


In eqn (4.75), the form of the checkerboard could be expressed more generally 
as 

pI = ki(-l)“ + fca (4.100) 

where is an arbitrary constant. A constant pressure field satisfying eqn (4.59) 
is called an hydrostatic pressure mode. For any contained fiow, the pressure is 
known up to an arbitrary constant and the hydrostatic pressure mode exists. The 
usual cure to that problem is to specify the pressure at one node. However in 
that case the constants ki and ki adjust so that the constraint is satisfied. But 
the checkerboard is not removed. In the treatment of the pressure mode, the 
averaging step only removes the nonconstant part of the checkerboard, and the 
final smoothed pressure is known up to a constant. 

It has been pointed out [18,19] that a method of the type described previously 
does not give good results at the boundaries and/or corners of the domain. More 
sophisticated algorithm including extrapolation procedures or special treatments 
at the boundary and corner nodes are described in [18,19]. 

4.5.5 The nine-node Lagrange element 

This element is standard for the mixed finite-element method using a quadra- 
tic velocity /bilinear pressure interpolation. It has been extensively used and its 
good properties have been long established [6,17]. 

This element contains 9 velocity nodes and 4 pressure nodes (Figure 4.1). The 
shape functions in the parent domain are given below: 

For velocity. 


= j(i - 0(1 - lUn 

(4.101) 


(4.102) 

VVa = -i(l + f)(l-D)£n 

(4.103) 


(4.104) 

JV, = j(l + £)(! + >()£>( 

(4.105) 

Afe = j(l -£’)(! + '»)') 

(4.106) 

JV7 = -j(1-£)(1 + x)£-; 

(4.107) 

JV. = -5(1 -£)(!->»’)£ 

(4.108) 

W. = (1 - £’)(! - 

(4.109) 
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The interpolation functions for the pressure are the same as for the four-node 
-H bubble element (see eqn (4.65)-(4.68)). 
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CHAPTER 5: NUMERICAL EXAMPLES 


5.1 Lid-driven cavity flow 

The lid-driven cavity flow is a typical example of internal flow. The simplicity 
of the geometry makes it one of the favored test problems in computational fluid 
dynamics. 

The motion of the upper wall of the cavity induces the development of the 
flow fleld. The flow being contained, viscous stresses are solely responsible for the 
initial motion of the fluid. However for high Reynolds numbers, the convective 
effects become significant and dominant for the development of the flow patterns,in 
particular for the occurence of recirculation regions at the bottom corners of the 
cavity. 

The Reynolds number is defined as 


Re = 


Uu,D 

u 


(5.1) 


where U^, is the tangential velocity of the upper wail, D is the size of the square 
cavity 2 uid u is the kinematic viscosity of the fluid. 

The domain and boundary conditions for both the primitive variable and the 
vorticity stream-function formulations are given in Figure 5.1. For the velocity- 
pressure formulation, the four-node + bubble element and the nine-node Lagrange 
element have been tested. For the later case as well as for the vorticity stream- 
function formulation, the meshes have 32 x 32 square elements, whereas a mesh 
of 20 X 20 square elements has been used for the four-node + bubble element. 
Moreover in the case of the four-node + bubble element, a Galerkin formulation 
was employed. 

For the velocity-pressure formulation, the boundairies are of pure Dirichlet 
type and are easy to implement. A pressure data is given at a point on the 
boundary to set the hydrostatic pressure level. For the vorticity stream-function 
formulation there are no boundary conditions for vorticity (because all the bound- 
aries are walls). The stream function is constant all over the boundeu'ies, and its 
normal derivative is zero for the left, bottom and right walls and is equal to the lid 
velocity at the upper boundary. In the vorticity stream-function formulation as 
well as the velocity- pressure formulation, the boundary conditions at the upper 
corners are discontinuous. However, this has produced no problem throughout the 
computations. 

The Reynolds numbers investigated in this study are 400 (Figure 5.2) and 
1000 (Figure 5.3) for the vorticity stream-function formulation and 400 (Figures 
5.4 and 5.5) for the velocity-pressure formulation. 

For this range of Reynolds numbers , the main structures of the flow are the 
primary vortex and the two secondary vortices. The primary vortex is created in 
the upper right region. As it grows, its center moves down toward the center of 
the cavity. The two secondary vortices develop at the lower left and right corners. 
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« wO,v"> 0 


b 


Figure 5.1. Lid-driven cavity flow : Problem domain and boundary conditions for 
the vorticity stream-function formulation (a) and the velocity-pressure formulation 
(b). 
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Figure 5.2. Lid-driven cavity flow: Steady-state solutions obtained by the vorticity 
stream-function formulation at Re=400; a) Streamlines; b) Corner streamlines; c) 
Vorticity. 
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Figure 5.3. Lid-driven cavity flow: Steady-state solutions obtained by the vorticity 
stream-function formulation at Re=1000; a) Streamlines; b) Corner streamlines; 
c) Vorticity. 
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b 


Figure 5.4. Lid-driven cavity flow: Steady-state solutions obtained by the velocity- 
pressure formulation using the nine-node element, at Re=400; a) Velocity; b) 
Pressure. 
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Figure 5.5. Lid-driven cavity flow: Steady-state solutions obtained by the veloc- 
ity pressure formulation using the four-node + bubble element, at Re=400; a) 
Velocity; b) Nonsmoothened pressure; c) Smoothened pressure. 
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For moderate Reynolds numbers, the span of the vortex at the lower right corner 
of the cavity is larger than that at the lower left comer. 

The results obtained for both formulations are in good agreement with avail- 
able data [19]. For a Reynolds number of 400, in both cases, the coordinates of 
the center of the primary vortex agree with the data given in [19] with a precision 
better than the node-to-node distance of the mesh. 

In the case of the four-node + bubble element, the pathology described in 
Section 5.4.3 , i.e, a pressure field polluted by a checkerboard mode, has been 
observed. The polluted pressure field has been treated through the procedure 
advocated in Section 5.4.4 to recover an acceptable pressure distribution. The 
results are shown in Figure 5.5. 

5.2 Plane jet flow 

This case has been studied only with the vorticity stream-function formula- 
tion. The problem domain and the boundary conditions are depicted in Figure 
5.6. The finite-element mesh has 64 x 64 square elements. For that problem, the 
Reynolds number is defined as 


Rc = 


V 


(5.2) 


where a is the width of the jet aperture, tJ is the mean velocity at the inlet and 
V is the kinematic viscosity of the fluid. The inlet velocity profile is parabolic and 
the Reynolds number is 167. The motivation for that simulation was to study 
the development of a plane jet in an open domain and our goal was to show the 
capability of our numerical scheme to simulate the breaking of symmetry in an 
initially symmetrical flow. 

In the initial state, the fluid is at rest. When the jet develops, the shearing 
effect between the front of fast moving fluid and the surrounding resting fluid 
induces the formation of two symmetrical vortices (Figure 5.7). The vortices grow 
as they move toward the right (downstream) boundary of the domain (Figures 5.8 
and 5.9). 

The selection of the boundary conditions for that case is important. The 
condition. 


dn 


= 0 


(5.3) 


which is applied on upper, lower and right (downstream) boundaries, implies 


Uf = 0 


(5.4) 


The effect of these boundary conditions on the flow field is especially seen in Figure 
5.9, where the streamlines at the downstream boundary curve and cross the limit 
of the domain at right angle. 
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Figure 5.6. Plane jet flow: Problem domain and boundary conditions for the 
vorticity stream-function formulation; a) Whole domain; b) Close view of the jet 
aperture. 
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Figure 5.7. Plane jet flow: Solution obtained by the vorticity stream-function 
formulation at t=1.2s; a) Streamlines; b) Vorticity. 
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Figure 5.8. Plane jet flow: Solution obtained by the vorticity stream-function 
formulation at t=3.6s; a) Streamlines; b) Vorticity. 
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Figure 5.9. Plane jet flow: Solution obtained by the vorticity stream-function 
formulation at t=4.8s; a) Streamlines; b) Vorticity. 
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After the two vortices associated with the eaorly development of the flow have 
left the computational domain, a steady-state synunetrical flow is obtained (Figure 
5.10). 

This flow is perturbed by modifying the inlet velocity profile over one time 
step. This is obtained by adding a 3% perturbation to the nodal values of the 
vorticity in the lower half of the inlet. The perturbed flow develops then in a 
nonsymmetrical pattern. Vortices form and pair as they are carried downstream 
(Figures 5.11 to 5.13). 

The amount of perturbation added to disturb the flow should not influence 
the final fiow field. Bristeau et al. [l] have reported results using a Galerkin for- 
mulation with a very fine unsymmetrical mesh. Their solution shows the breaking 
of symmetry without adding a perturbation. The use of a nonsymmetrical mesh 
acts as a perturbation that triggers the development of a nonsymmetrical solution. 

5.3 Flow past a circular cylinder 

The fiow of a viscous fluid past a circular cylinder is one of the most extensively 
studied cases and has become a classical test for numerical procedures [5,8]. It 
is of engineering interest because it presents most of the flow features that are 
encountered in more complicated external flow situations such as separation and 
vortex shedding [20,21]. In fact the cylinder can be seen as the basic shape that 
can be degenerated in all the more complicated symmetrical and nonsymmetrical 
profiles. This case has been investigated using both the primitive variable and the 
vorticity stream-function formulations. In both cases, symmetrical flow solutions 
are first obtained and then perturbed to produce vortex shedding in the wake of 
the cylinder. 

The meshes and boundary conditions for the two formulations are shown in 
Figures 5.14 and 5.15. For the velocity-pressure formulation, only the nine-node 
Lagrange elements were used. 

For this flow, the Reynolds number is defined as 

Re = (5.5) 

1 / 

where Uoo is the inlet velocity, d is the diameter of the cylinder, and i/ is the 
kinematic viscosity of the fluid. The present study has been carried out for Re = 
100 * 

The specified inlet velocity profile is a constant. For the upper and lower 
boundaries, in the velocity-pressure formulation, we specify zero tangential stress 
and zero vertical velocity whereas in the vorticity stresmi-function formulation, the 
value of the stream function is specified and the vorticity is set to zero. From a 
physical point of view, these boundary conditions are equivalent. However they are 
only valid if the boundaries are far enough from the cylinder so that the vorticity 
generated in the boundary layer at the cylinder surface does not diffuse up to the 
outer boundaries. In a finite-difference framework, the boundaries of the domain 
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Figure 5.11. Plane jet flow: Solution obtained by the vorticity stream-function 
formulation at t=1.2s after perturbation; a) Streamlines; b) Vorticity. 
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Figure 5.12. Plane jet flow: Solution obtained by the vorticity stream-function 
formulation at t=2.4s after perturbation; a) Streamlines; b) Vorticity. 
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b 

Figure 5.13. Plane jet flow: Solution obtained by the vorticity stream-function 
formulation at t=3.6s after perturbation; a) Streamlines; b) Vorticity. 
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a 





b 


Figure 5.14. Flow past a circular cylinder: Problem domain and boundary con- 
ditions for the vorticity stream-function formulation (a) and the velocity-pressure 
formulation (b). 
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Figure 5.15. Flow past a circular cylinder: Finite-element meshes for the vorticity 
stream-function formulation (a) and the velocity-pressure formulation (b). 


46 





can be set at an arbitrary large distance from the cylinder through mapping tech- 
niques. In the finite-element method, procedures such as exponential stretching 
are used to render the infinity of the computational domain. In the present study, 
following the work of Gresho [5] and Hughes [7], we used a domain of finite extent. 
This presumably induced a blockage effect in the fiow around the cylinder [7], but 
did not prevent us from simulating the main fiow features such as the recirculation 
and the vortex shedding. At the downstream boundary, the stress is set to zero in 
both tangential and normal directions in the primitive variable formulation and 
the normal derivatives of vorticity and stream function are set to zero in the vortic- 
ity stream-function formulation. Those boundary conditions are not equivalent. 
Setting the stress to zero imposes no constraint on the direction of the fiow at 
the boundary whereas a zero normal derivative of the stream function forces the 
streamlines to leave the domain in a direction normal to the boundary. However, 
this effect is of small influence on the computed fiow. At the cylinder surface, the 
boundary conditions are naturally the no-slip conditions which are staightforward 
in the velocity-pressure formulation. For the vorticity stream-function equations, 
as explained earlier, there are no obvious boundziry conditions at a solid surface. 
In the present case, we take advantage of the symmetry of the domain to set the 
value of the stream function at the cylinder. When the nonsymmetrical fiow is 
studied, this condition is no more valid. In fact because the stream function is 
set at the upper and lower boundaries and at the cylinder surface, the fiow rate is 
constrained to be the same above and under the cylinder. The previous remarks 
show that for external flows, the boundary conditions of the primitive variable 
formulation are more accurate and easier to implement than those of the vorticity 
stream-function formulation. 

For both cases, the flow develops in a symmetrical pattern (Figures 5.16 and 
5.17). Two standing vortices grow in the wake of the cylinder. In experimental 
studies [22] of such a flow, nonsynunetrical solutions presenting vortex shedding 
have been observed for Reynolds number greater than 40. Symmetrical solutions 
can not be obtained experimentally because the flow is unstable for any pertur- 
bation, and of course, in an experimental set up perturbations cannot be avoided. 
This is not the case for a numerical study. The mesh we have used is symmetrical 
and all the conditions are symmetrical up to the last digits. The computations 
being performed in double precision, the level of unsymmetry is virtually reduced 
to zero. 

The nonsymmetrical flow is obtained by perturbing the steady symmetrical 
solution. In the case of the vorticity stream-function formulation, this is done 
by introducing some artificial vorticity at the centerline in the near wake of the 
cylinder before restarting the computations. The perturbation amounts to 5% of 
the maximum value of the vorticity in the flow field. For the primitive variable 
formulation, the cylinder is rotated in the counterclockwise direction over four 
time steps. It would be simpler to perturb the velocity field at the centerline in 
the wake of the cylinder, but this would produce a nonsolenoidal flow field which 
could lead to nonconverging computations. 
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c 


Figure 5.16. Flow past a circular cylinder: Steady-state solutions obtained by 
the vorticity stream-function formulation at Re=100; a) Streamlines; b) Stationary 
streamlines; c) Vorticity. 


b 


Figure 5.17. Flow past a circular cylinder: Steady-state solutions obtained by the 
velocity-pressure formulation; a) Velocity; b) Pressure. 
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In both situations the perturbation produces a nonsymmetrical flow field. The 
symmetry of the standing vortices is broken. They start oscillating and finally shed 
(Figures 5.18 to 5.23). 
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Figure 5.19. Flow past a circular cylinder: Solutions obtained by the vor- 
ticity stream-function formulation at Re=100 at t=240s after perturbation; a) 
Streamlines; b) Stationary streamlines; c) Vorticity. 
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Figure 5.20. Flow past a circular cylinder: Solutions obtained by the vor- 
ticity stream-function formulation at Re=100 at t=360s after perturbation; a) 
Streamlines; b) Stationary streamlines; c) Vorticity. 
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Figure 5.21. Flow past a circular cylinder; Solutions obtained by the velocity- 
pressure formulation at Re=100 at t=15s after perturbation; a) Velocity; b) Pres- 
sure. 
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Figure 5.22. Flow past a circular cylinder: Solutions obtained by the velocity- 
pressure formulation at Re=100 at t=105s after perturbation; a) Velocity; b) 
Pressure. 
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b 


Figure 5.23. Flow past a circular cylinder: Solutions obtained by the velocity- 
pressure formulation at Re=100 at t=180s after perturbation; a) Velocity; b) 
Pressure. 


56 





CHAPTER 6: CONCLUSIONS 


The present study has focused on the development of streamline-upwind/ 
Petrov-Galerkin finite-element procedures for the Navier-Stokes equations in both 
the primitive variable and the vorticity stream-function formulations. The abilities 
of the two methods for simulating internal and external fiows have been studied. 
For all test cases the results are in good agreement with available data. 

The main advantages of the vorticity stream-function formulation are the 
simple form of the equations in two-dimensions and the in-built satisfaction of the 
incompressibility constraint. 

In two dimensions, the vorticity transport equation and the Poisson equation 
for the stream function are scalar and there are only two degree of freedom in 
the problem. Moreover, the vorticity stream-function form of the Navier-Stokes 
equations allows equal order of interpolation for the vorticity and the stream func- 
tion. In fact the bilinear interpolations which have been used here for both the 
unknowns are sufficient which is an asset from the point of view of implementation. 

On the other hand, a fiow field obtained by the solution of the vorticity 
stream-function equations is by definition divergence-free and the initial condition 
of the problem do not need to satisfy the incompressibility constraint. This allows 
the use of an initial fiow field noncontinuous at the boundary of the domain. 

In the case of the primitive variable formulation in two dimensions, the num- 
ber of degrees of freedom is three. Moreover in this formulation the use of equal 
orders of interpolation for pressure and velocity can lead to erroneous pressure 
solutions such as checkerboard modes. Mixed interpolations are usually used in 
which implementation is less straightforward than equal order of interpolation. 
Here we implemented a mixed interpolation using the nine-node Lagrange ele- 
ment and a four-node + bubble element. The latter exhibits a checkerboad mode 
which can be treated through smoothing techniques. On the other hand, in the 
velocity-pressure method, the pressure is adjusted to satisfy the incompressibility 
constraint in a weak sense. 

From the point of view of the boundary conditions, the velocity-pressure for- 
mulation allows a more natural implementation than the vorticity stream-function 
formulation especially for fiow domains having inner boudaries. In that case,the 
vorticity stream-function formulation provides only approximate boundary condi- 
tions for the stream function and no conditions for vorticity. In the present study, 
the simulation of the fiow around a cylinder is possible under the assumption of 
symmetry of the fiow pattern at the upstream of the cylinder. When the fiow is 
perturbed, this assumption is not valid, but in the present scheme, the streaun 
function cannot be updated at the inner boundary. The case of a fiow with mul- 
tiple inner boundaries, such as an array of cylinder, cannot be handled by this 
scheme. 

The velocity-pressure formulation allows a good treatment of the boundary 
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conditions. The no-slip and no-penetration conditions at solid surfaces are easy 
to implement. Moreover there would be no problem for handling flows in domains 
with multiple solid boundaries. 

At a last point, the extension of the codes to three dimensions is not easy 
for the vorticity stream-function formulation. In three dimensions, the vorticity- 
stretching term does not vanish and the equation system has two nonlinear terms. 

The primitive formulation can be extended to three dimensions without major 
difficulties. In the present study, the use of an explicit time-integration scheme 
obtained via a mass-lumping technique reduces the storage requirements of the 
scheme, and would be a significant asset for the solution of three-dimensional 
problems. 

The possibilities for improvements and further developments of the present 
work are manifold. 

In the case of the vorticity stream-function formulation the main improve- 
ment that can be made is a better treatment of the boundary conditions for the 
vorticity and the stream function at solid surfaces. In particular the need for a 
scheme allowing the treatment of flows in domains with multiple solid boundaries 
is obvious. This has been implemented by Liou [24] following the scheme proposed 
in [1]. 

Two major directions of development are possible for the velocity-pressure 
formulation. 

The first possibility is to add the temperature as a physical unknowns in the 
problem. This can be obtained by coupling the energy equation to the momentum 
equation and the incompressibility constraint. Based on the present formulation, 
Ganjoo [25] has implemented that problem and obtained results comparing well 
with those presented in [5|. 

The second major development to the present work is the extension to three 
dimensions. As described earlier, the velocity-pressure formulation has been im- 
plemented in order that this extension be possible without dramatic changes in 
the algorithm. The trend in CFD is to go to three-dimensional simulations, and 
the present formulation is a good basis for such development. 
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